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Forward 


The  first  two  parts  of  this  report  wind  up  a  few  questions  in 
the  mathematical  formulation  of  vector  fields  governed  by  conser¬ 
vation  and  rotationality  laws,  with  explicit  application  to  fluidyn- 
amic  fields,  possibly  with  shock  waves.  The  points  treated  have 
a  strong  bearing  on  computational  schemes  and  the  stability  of 
numerical  calculations  and  the  results  provide  a-priori  informa¬ 
tion  on  the  way  to  select  the  appropriate  set  of  equations,  the 
right  functional  and  the  most  promising  approximation  space  for  a 
finite  element  descretizations .  The  last  assertion  is  then  tested 
for  the  tricomi  equation  in  a  non-uniformly  elliptic  domain. 

A  mixed  Tricomi  problem  is  descretized  by  an  alternative 
collocation  scheme  which  proves  to  be  accurate  and  stable  as  dem¬ 
onstrated  on  a  few  test-cases.  The  collocation  finite  difference 
schemes  have  superior  to  the  finite  elements  ones  for  this  case, 
so  far.  They  are,  however,  more  specialized  and  natural  for  linear 
problems  and  simple  geometries.  The  variationally  based  finite 
elements,  on  the  other  hand,  hold  a  better  promise  for  complex 
geometries,  and  for  an  accurate  treatment  of  shocks. 


I.  A  Note  on  how  to  select  from  too  many  equations  the  right  ones 

to  solve  a  given  problem. 


Nima  Geffen 


Abstract 

The  field  laws  for  a  physical  continuum  are  often  described 
by  too  many  first  order  partial  differential  equations  for  the 
number  of  required  field  quantities.  The  question  described  and 
a  simple  way  to  resolve  it  is  given  for  the  so-called  conserva¬ 
tive  and  non-conservative  representations  of  continuum  mechanics 


General  Equations 

Continuum  mechanics,  electrodynamics  and  other  physical 
theories  can  be  modelled  by  various  specialization  of  the  follow¬ 
ing  equations: 

Cl)  VQ  A  =  G 

C2)  v  x  u  =  W  Cn  equations) 
for: 

x  =  x^  i  =  l,...,m  independent  variables 
u  =  Uj  j  =  l,...,n  dependent  variables 
A  Cx^,Uj)  =  A^  k  =  l,...,m 

and: 

GCxi,u.),  W(xi,uj)  =  W<j) 

where  the  source  function  G  is  arbitrary,  but  the  vorticity  W 
has  to  satisfy  a  compatibility  condition: 

(3)  V0  W  =  0 

The  system  (1)  (2)  includes  (n+1)  first-order  partial 
differential  equations  for  the  n  unknown  u  j .  The  overdetermin- 
acy  is  apparent  only,  due  to  the  fact  that  the  n  rotationality 
conditions  are  not  independent  (note  eq.  (  3)  )  because  any  (n-1) 
statements  imply  the  one  left  as  will  be  shown  explicitly  in  the 
following . 
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Examples 

1.  Maxwell’s  equations  for  the  electromagnetic  field  are  [  :j 

el)  V-H  =  0  e3)  VxE  =  -  -  H. 

“  —  C  — ' t 

e2 )  V  •  E  =  4irp  e4 )  VxH  ~  —  E  +  —  2 

—  —  c  — t  c 

Equations  el)  -  e4)  are  not  independent:  for  smooth  xieids 
(twice  differentiable)  eq.  3.  implies  (V-H)^  =  0  and  eq.  el) 
holds  for  all  times  if  it  holds  for  t  =  tQ,  thus  it  is  an  ir.it  i 
condition  (at  most). 

The  electric  charge  density  and  currents  (p  and  j  respect 
ively)  cannot  be  prescribed  arbitrarily,  since  eq .  e2  and  e*t) 
imply  a  constraint  on  their  source  terms  p  and 

e5)  pt  +  |  V-j  =  0  (continuity  of  charge) 
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II.  Steady-state  Aerodynamics 
for  which: 

x  =  (x^,X2»Xg)  =  (x,y,z)  -  space  coordinates 
u  =  (u^jUjjUg)  ~  Cu,v,w)  -  velocity  components 
A  =  p  u 


and : 


p  =  p(u  )  density 


G  =  0,  W  =  W( x ) 

where:  W  -  0  for  irrotational  flow,  but  changes  across  o.n'e 

shocks.  Spelled  out  in  Cartesian  coordinates  we  get: 


2  2  2  2  2  2 

e6)  (a  -u  )u  +  (a  -v  )v  +  (a  -w  )w  -  2uww  -  2uv  V  —  2  7WW 
x  y  z  x  y  2 


( i)  w  -  u  =  W 
x  z 


(1) 


(ii)  w  -  v  =  W 

y  z 


(2) 


(iii)  Uy  -  vx  =  W 


(3) 


i.e.  4  equations  for  the  3  unknown  functions  Cu,v,w). 
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Methods  of  Solution  and  Descretization 

The  system  equations  (1),  (2)  even  when  simplified  for 
specific  physical  field,  is  multidimensional  and  coupled,  which 
renders  it  complicated  to  analyse  and  inconvenient  to  solve. 
Auxiliary  functions  (e.g.  scalar  and  vector  potentials) 
have  been  taylored  to  simplify  and  clarify  the  mathematical 
picture  and  render  it  elegant,  intelligible  and  solvable. 

For  irrotational  fields  u,  i.e.  W  a  0  eq.  (1)  reduce  to 
one  second  order  equation  for  the  scalar  potential  <f> ,  defined  by 
u  =  V4>,  so  that  eq.  (2)  is  satisfied  identically. 

For  Maxwells  equations  the  standard  analysis  is  done  via 
the  scalar  and  vector  potentials:  <j>  and  A  respectively,  (e.g.  L 
where : 

E  =  ‘  c  It  *  '  Vx-' 

Written  in  terms  of  (<f>,A)  Maxwell’s  system  reduce  to  4  equations 
for  the  4  components. 

The  resulting  equations  are  higher  order,  and  admit  a  wider 
family  of  solutions,  all  equivalent  under  gauge  transformations: 

ft(x,t) 

A'  =  A  +  Vf 
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or  in  4  dimensional  notation: 

(A =  (A,-*)  +  [Vf,-  |  ft3 

The  potential  formulation  for  the  electromagnetic  field  is 
endowed  with  a  beautiful  structure  ,  lucid  transformation  prop¬ 
erties  and  striking  accessible  information  content  Ce.g.  de¬ 
coupling  into  2—  order  wave  equations  for  each  of  the  components) 
unfolding  the  wealth  of  electromagnetic  waves  ref.  [1]. 

Another  suggestion,  to  equalise  the  number  of  equations  and 
unknowns  is  to  add  a  new  dependent  variable: 

v  =  V  •  u_ 

and  using  the  relation: 

2 

VxVxu  =  V ( V • u)  -  V  u 

replace  eq.  ( 2 )  by  its  rotor: 

(4 )  72u  =  Vv  -  VxW, 

which,  with  eq.  (1)  gives  (n+1)  equations  for  the  (n+1)  unknowns 
(u , v) . 

The  formulation  above  has  been  suggested  by  M.  Mock  for  compu¬ 
tational  purposes,  with  a  staggered  mesh  for  (u,v)  Cto  avoid  decoup¬ 
ling  of  the  descretized  equations  for  u  and  v)  to  solve  boundary 
value  problems. 
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Direct  field  formulation 

The  auxiliary  function  is  (e.g.  potentials)  formulations 
invariably  raise  the  order  of  the  equations  to  be  solved;  the 
first-order  system  becomes  second-order.  This  requires  a  higher 
degree  of  smoothness  for  the  solution  function  and  may  be  a 
draw-back  for  numerical  analysis  and  calculations.  Thus,  although 
many  large  computer  simulations  are  based  on  'potential'  formu¬ 
lations  (e.g.  ref  [2])  a  direct  solution  of  the  first  order  sys¬ 
tem  has  been  found  beneficial  f  3 ] .  Sometimes  essential  [4], 
-especially  for  'initial'  rather  than  boundary  value  problems. 

The  question  is  how  to  choose  the  'right'  (n-1)  rotationality 
conditions  that  will  give  with  the  continuity  eq.  (1),  the  'right' 
(nxn)  system  for  a  stable  descretization  for  a  marching  scheme 
to  solve  the  first  order  system,  to  obtain  directly  the  n  field 
components  (so-called  primitive  variables).  A  simple  treatmen:. 
for  irrotational ,  3-dimensional  fields  (worked  out  for  the  prob¬ 
lem  in  [4])  is  given  in  [  5] .  It  is  extended  here  for  the  generd. 
case  teq.  (1),  (2)),  where  the  field  may  have  sources  and  be 
rotational  (e.g.  an  electromagnetic  field  with  moving  charges, 
flow  behind  a  curved  shock,  motion  of  reacting  gases.) 
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Choice  of  Equations  for  3-Dimensions 

Spelled-out  for  3  dependent  variables  and  3  independent 
ones  (e.g.  Cartesian  space  coordinates)  we  get: 


x  =  (x^,x0 ,x J  =  (x,y,z) 
u  =  (u-L,u2>u3)  =  (u,v,w)  (x,y,z) 

A  =  (A(1) ,A(2) ,A(3))  (u,v,w;  x,y,z) 


and  the  system  of  equations  is: 


Cl)  AC1)  +  A(2)  +  A(3)  =  G 
x  y  z 


w  -  v  =  W 

y  2 


(i) 


(  2  )  u  -  w  =  W 
z  x 


,(2) 


v  -  u  =  W 
x  y 


(3) 


i.e.  4  equations  for  the  three  unkonwon  functions  (u,v,w),  a 
redundant  system  for  a  well  defined  field. 

In  addition,  a  compatibility  on  W  requires: 


C3) 


W(1)  +  W(2)  +  W(3)  =  0 
x  y  z 
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Irrotational  Fields 

In  this  case  W  =  0,  eq.  (3)  is  satisfied  automatically  and 
equations  (2)  become: 

i)  w  -  v  =  0 

y  z 

ii)  u  -  w  =0 

z  x 


iii)  v  -  u  =  0 
x  y 


for  (u,v,w)  (x,y,z)  unique  and  sufficiently  smooth,  (e.g. 
twice  differentiable  in  each  of  the  independent  variables)  we 
get : 


ii)  —  ■>  0  =  u2y  -  wxy  =  Cuy)z  -  (wy>x  =  (vx>z  -  (wy)x  =  <v2-wy)x 


9 / 3x  (ii)  (u,w)  €  C‘ 


(i) 


ve  ( 


,2 


v  -  w  =  c(y  ,z) 

z  y 

For  u  irrotational  at  any  x  =  xn , 


(vz-Wy)  (xQ  ,y  ,z)  =  0  •*  C(y,z)  =  0  c(y,z) 


and  i),  ii)  >  (iii)  which  can  be  considered  and  "initial  condition" 
at  most,  (e.g.  x^  — >  ,  and  a  uniform  field  there).  In  exactly 


the  same  manner: 
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i)  iii)  — >  [u  -w  ]  =  0 

z  x  y 

u  — w  =  F(x,z) 
z  y  ’ 

and  irrotationality  at  any  y  =  y^  plane  implies  F  =  0  and  eq.  (ii) 
and: 

i )  i  i )  — >  [  v  -  u  J  =0 
x  y  2 

and  irrotationality  at  any  z  =  zQ  implies  iii). 

Thus  any  of  the  3  irrotationality  conditions  implies  the 
third,  for  unique  smooth,  irrotational  fields  everywhere.  For 
stable  numerical  algorithms  however,  one  has  to  choose  set  that 
carries  the  appropriate  "boundary  information"  along  the  march¬ 
ing  coordinate;  if  the  integration  is  to  be  carried  out  along 

-A_ 

the  x^  direction,  the  i—  component  of  the  iirrotationality 
equations  has  to  be  omitted  [4]  ,  [5], 
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For  the  general  irrotational  Case; 

*  .  .  .  „(1) 

(1)  w  -v  -W 

y  2 


Cii)  az  “  wx  '  w 


Ciii)  vx  "  uy  =  W 


(3)  W(1)  +  W(2)  +  W*3)  =  0. 

x  y  z 


Differentiating  and  substituting  we  get 


-2.  i)  w  -  v  ^ 

3x  yx  x  x 


(w  )  -  v 

x  y  zx 


Cu)  -  W(2)  -  v  =W(] 
z  y  y  zx  x 


Cu  -v  )  =  =  _W<J» 

y  x  z  x  y 


(  iii) 


v  -  u  =  W(  3)  +  F(x,y) 
x  y 


[(vx-uy)  -  wc 3)](x,y,zQ)  =  F(x>y»zo) 

Thus  i)  ii)  — iii)  provided  vx~uy  is  given  at  z  =  zQ . 

In  the  same  manner  each  two  of  the  equations  determine  the 
third,  for  which  initial  conditions  have  to  be  given  for  uniquene 
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Algebraic  treatment 

Equations  Cl),  (2)  can  be  written  in  a  matrix  form  (used 
widely  for  numerical  analysis  and  descretizations ) : 


(5 ) 


with  the  corresponding  (4x3)  matrices  C  and  4  forcing  functions  F. 


k 


-  13  - 

The  system  (.51  is  again  redundant,  and  one  of  the  last  3 
equations  has  to  be  deleted  to  render  the  (4x3)  coefficient 
matrices  C  into  (3x3)  matrices  £.  The  choice  is  directed 
by  the  marching  "time  like"  direction,  for  which  the 
matrix  has  to  be  invertible,  hence  nonsingular,  hence  with 
no  row  (or  column)  of  zeros.  This  automatically  rules  out  one 
equation:  when  (u,v,w)  are  given  at  x  =  Xq  (e.g.  xQ  — > 
for  steady  flow  about  an  obstacle),  the  first  irrotationality 
condition  has  tc  omitted.  A  marching  procedure  along  the  x 
axis  requires  the  inversion  of  the  non-redundant  matrix  £(x) 
such  that: 

uv  --  -[fU)]-1  i'y’u,  >  ♦  c61x,]_1f  < 

x  y  Z 

In  the  same  manner,  integration  schemes  along  the  y  and  2 
directions  require  the  omission  of  (5i0  and  (5iii),  respectively, 
as  a  necessary  condition  (not  always  sufficient!)  for  a  well 
defined,  stable  scheme,  (as  is  seen  in  [ 4 J ) . 
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Final  Remarks 

The  question  and  system  treated  are  elementary  and  so  gen¬ 
eral,  the  analysis  and  answer  so  simple, to  be  judged  trivial, 
if  not  for  the  fact  that  the  question  does  come  up  occasionally 
Ce.g.  [4]),  and  the  answer  not  always  immediate.  Essentially 
the  same  problem  has  been  treated  recently  (and  has  come  to 
our  attention  while  writing  this  note),  in  a  completely  differ¬ 
ent  context  for  different  reasons  and  aims  in  [63.  It  also  seems 
to  be  related  to  other  recent  approaches  [ 7]  and  to  variational 
formulation  and  analysis  of  the  system  at  hand  with  related 
physical  applications  [8]. 

Work  on  relaxing  the  smoothness  requirement,  re-formulation 
and  analysis  for  non-smooth  fields  (e.g,  flows  with  shocks)  is  in 


progress . 
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II.  ALTERNATIVE  VARIATIONAL  FORMULATIONS 


FOR  NONLINEAR  VECTOR  SYSTEMS 


Nima  Geffen 


Abstract 

Variational  formulations  for  vector  fields  described  by  a 
system  of  partial  differential  equations  of  whatever  type,  possib¬ 
ly  nonlinear,  and  with  initial-boundary  conditions  are  given. 
Smoothness  properties  of  suitable  approximation  spaces  are  viewed 
and  the  effect  of  coercive  ,  natural  boundary  conditions  dis¬ 
cussed  briefly.  Examples  are  drawn  from  electromagnetic  field 
theory  and  f luidynamics . 


Differential  equations 

Consider  the  following  conservation  and  rotationality  condi¬ 
tions,  for  a  vector  field  u(x) : 

(1)  V-V  =  G 

(  2  )  V  x  u  =  W 

(з)  uCx)  =  f (x)  x  £  3fi: 
where : 

u  =  _u,.  is  the  vector  of  dependent  variables  i  =  l,...,::* 

x  =  Xj  is  the  vector  of  independent  variables  j  =  l,...,n 

G  =  G(x)  is  a  given  function 

W  =  Wj (x)  is  the  vorticity 

V  -  V^(x,u) 

The  vorticity  W  has  to  obey  a  compatibility  condition: 

(и)  7 • W  -  0 

The  system  of  equations  (1),  (2)  is  quite  general,  it  can 
be  linear  and  nonlinear,  elliptic,  hyperbolic  or  mixed  with  smooth 
or  non  smooth  solutions.  The  independent  variables  x^  may 
designate  space  and  time  coordinates  and  different  kinds  of 
initial  and/or  boundary  conditions  may  be  appropriate  for  differ¬ 
ent  problems.  Higher  order  equations  may  be  put  into  this  form, 
and  examples  of  applications  include  the  description  of  electro¬ 
magnetic  field  s,  theory  of  elasticity,  f luidynamics ,  and  plasm, 
dynamics,  including  flows  with  shocks. 


2 


Variational  Formulations 

A  variational  formulation  of  the  field  u  satisfying  (1), 
(2),  (3)  is  a  functional  J(v)  defined  on  Q  whose  stationary 
value  is  obtained  for  v  =  u: 

<5J(v)  =  0  <*=v=u 


For  well  posed  problems,  for  which  u(x)  is  unique 


6J(v)  =  0  «—*  v  =  u 


Variational  formulations  are  scalar,  short,  additive  (the 
functionals  for  complex  systems  are  direct  sums  of  their 
simpler  parts),  invariant  under  appropriate  classes  of  trans¬ 
formations  and  are  often  convenient  for  theoretical  analysis 
and  for  numerical  simulations,  e.g.  by  the  finite  elements 
method.  Integrals  can  easily  be  descritized  and  approximated, 
and  the  smoothness  requirements  on  the  functions  are  less  strin¬ 
gent  than  for  the  corresponding  differential  system.  This  last 
point  is  most  important  from  the  numerical  view  point  in  addition 
to  a  better  rationale  for  the  treatment  of  shocks. 

The  case  G  =  0  is  described  in  [1].  The  functional  J(v,x): 


J 


f 

L(y,x)  +  X.CVxv-W)  + 

« 


Xxv>  do 


an 


m 


is  stationary: 
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<5J(v)  =  0  at  v  =  u 

provided  that:  V  x  v  =  0. 

u  — 

V  =  V  L 
—  u 

resulting  in: 

V  =  -  VxA_ 

The  variation  is  done  on  all  v  for  which  J  is  defined 
and  v  satisfying  coercive  boundary  conditions  on  the  boundar 
3ft^  =  or  for  which  A_  ||  v  or  v||do_  on  3ft. 

For  the  non-sourceless  (or  sourceful)  case;  the  foliowin 
variational  statements  hold: 

Th.  1 

The  functional: 

(5)  J(v)  =  J[L-g* 

n 

is  stationary  for 

(6)  V  x  V  -  0 

u  — 

where : 

(7)  V-g  =  G  or  £  =  V 

(8)  V  =  V  L 

—  u 


v  +  Ay  (  Vxv-\V)  ]  •  dx  +  j  Axv  »  dp 

m 


v  =  u  satisfying  Cl),  (2),  (3)  provide; 
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C9)  V  =  VxX_  +  g 

v  is  allowed  to  vary  over  all  functions  for  which  (5)  is 
defined  and  finite  and  which  satisfy  the  coercive  B.C.  (3). 
The  proof  is  straight  forward  and  follows  the  details  in  [1] 
exactly . 
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Corollary 
The  functional: 


(10)  J ( v) 


j  (L-£*  V  + 

ft 


X»W)dx  + 


Xxv*  do 


3ft. 


m 


(3) 


is  stationary  for  v  =  u  satisfying  (1),  (3)  when  the  variation 
is  taken  over  all  fields  satisfying  eq.  (2)  and  the  initial/bound¬ 
ary  conditions  (3):  i.e. 


6  J  (  v) 


0  <—  v  =  u  :  7-V  =  G 

yxV  =  W 


v.  =  £.  (x)  x£3n.- 

J  J 


The  restricted  variational  statement  (10)  follows  from  the 


statement  in  (5)  by  inspection 


An  alternative  formulation  is  obtained  by  integrating  the 
2  term  by  parts,  using  the  vector  identity: 

y-(Xxv)  -  v*  VxA_  -  A_*  Vxv 

substituting  in  (5): 


Th.  2 


X*Vxv  =  v.VxX 


The  functional 


V • ( Xxv) 


(11)  J(v)  =  j  (b-£*v  +  v-VxX_  -  X  •  M)dx 
is  stationary  for 


v  =  u  satisfying  (1),  (2). 

In  the  variational  formulation  (11)  is  made  over  all  v 
in  Q  which  render  all  terms  integrable.  The  lagrange  multi¬ 
plier  X_  is  required  to  have  integrable  first  derivatives,  which 
appear  explicitly  in  J. 

The  surface  term  drops  out,  and  the  solution  v  =  u  satis¬ 
fies  the  natural  boundary  conditions:  X.IIH  or  un  do 


or.  3  S' 


Smoothness  requirements 

In  the  variational  formulation  (5)  v  is  required  to  be 
at  least  once  differentiable  (for  the  2—  term  to  be  defined) 
which  holds  also  for  (10)  (where  eq.  (2)  has  to  be  satisfied). 

The  lagrange  multiplier  _A  in  this  formulation  can  be  just  inte- 
grable,  e.g.  a  step  function.  The  variational  statement  (11), 
on  the  other  hand,  does  not  involve  derivatives  of  v  (hence 
admit  integrability  -  only  there)  but  requires  £  to  differ¬ 
entiable  at  least  once.  Formulations  (5),  (10)  include  a  sur¬ 
face  term  and  involve  spaces  satisfying  appropriate  boundary 
conditions;  in  statement  (11)  the  surface  term  has  dropped  out, 
and  the  solution  v  =  u  making  it  stationary  satisfies  a  natural 
boundary  condition.  Both  smoothness  requirements  and  behavior  on 
the  boundary  has  bearings  on  approximation  spaces  used  for  numer¬ 
ical  calculations,  and  on  the  stability  of  numerical  schemes. 

Simple  examples  are  described  in  the  following  part  of  this  report. 
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Examples : 

1.  Steady  electromagnetic  field 

Maxwells'  equations  for  a  steady  state  can  be  written  av  [?] 
i)  V*E-p  ii)  VxE  -  0 

iii)  V*B_  -  0  iv)  VxB_  =  j 


A  variational  statement  for  i)  ii)  is: 


1  7  ( F\ 

( 5 '  )  JCE)  =  j (E  '2-£*E)  +  Ad  ' ( VxE)  + 


A(E)xE-do 


. 


for  E  arbitrary,  or: 


CIO')  J(E)  -  j  (EV2  -  g'E)  for  E  irrotaticnal . 

h 

where:  ^  is  any  solution  to 

=  p  or:  £  =  V-1p 

The  variational  statement  (11)  for  i)  ii)  become: 


(11')  JCE) 


[E2/2-£-E)  -  E-CVxft] 


The  corresponding  variational  formulations  for  the  magne 


equations  ixi)  iv)  are: 
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(5")  J(B)  =  f  B 2 /  2  +  \_(B)(VxB-j  )  +  f  /B)xB-da 
J  i 

dQn 

(10")  J  ( B  )  -  I  B 2  /  2  -  a_B  •  j  +  ji(B)xB-dc 

for  £  satisfying  eq ,  iv) 
and 

til”)  JCB)  =  j  L  E 2  /  2  -  B.(VxACB))] 

the  functionals  for  the  combined  field  are  obtained  by  simply 
adding  the  ones  for  the  'separated'  system: 


(5'")  J(E,E)  = 


rTE2-B2 


,  ,B  . 

2  -  £*L  +  i  ‘2 


+  _A(E)(VxE)  -  _A(B)(Vx3)]  + 


AExE • do 


an 


rr. 


r  u 

!  A  xB • do 


an 


(10'")  J (  B , E )  = 


E  2  -  B  2 


in 


-  £•  E  +  a_B  •  3 


an 


for  irrotational  E  and  B  satisfying  iv) . 

The  statement  (10")  can  be  reduced  to  the  one  used  for  the  scala 
and  vector  potentials  for  the  irrotational  elective  and  solenoid 
al  magnetic  fields,  (e.g.  [3]  pg.  366  eq.  (11-65). 
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Finally  (11)'  and  (11")  combine  to  give: 


(11"'  ): 

J(E,B) 


-  g*E  -  E* VxXE 


+  B-Vx_XB 


i 

i 

F 
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2.  Steady  fluidynamics 

The  differential  equations  are: 

V* (pq)  =  0  P  =  P (q^  ) 


Vx*<^  =  w 


and  the  corresponding  Lagrangian  L  and  functionals  are: 


ch 


or:  pu. 


J(u)  -  j  [  L  +  X.*(^xu  -  w)  ]  +  I  Xxu*do 
U  3  Q 

(lfr)  or: 

f  f 

JCu)  ^l^-X^w  +  j  Ax  u*  da 
0  3fi 

or: 


Cll)  ,J(u)  =  [  L  -  k* w  +  u*7xA_ 

Q 

for  the  appropriate  function  spaces,  with  the  required  smooth¬ 
ness  properties ,  and  constraints  in  the  region  and/or  on  the 
boundary . 

The  results  <$)  have  been  described  and  analyzed  in  [la],  [lb] 
where  applications  to  specific  flows  are  given;  the  formulation 
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(lft)  of  a  slight  modification  used  in  electromagnetics  and  other 
field  theories  but  not  generally  useful  for  computations  due  to 
the  practical  difficulty  in  actually  constructing  a  wide  enough 
family  of  v  fields  that  satisfy  the  rotationality  condition  (2)). 
The  formulation  (ll)  is  new  and  can  be  readily  specialized  to  par¬ 
ticular  fields;  it  offers  an  interesting  alternative  from  the 
theox',eti cal  point  of  view  but  holds  rather  doubtful  premise  compu¬ 
tationally  (stability  problems?). 
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Concluding  Remarks 

A  preliminary  report  on  alternative  variational  formulations 
is  given  for  general  vector  fields  governed  by  systems  of 
partial  differential  equations,  possibly  nonlinear  specifying 
their  sources  (eq.  Cl),  e.g.  conservation  of  mass  for  G  =  0) 
and  vorticity  Ceq.  (2),  e.g.  w  =  0  for  irrotational  fields). 

This  framework  is  pregnant  with  information  and  connections 
with  other  variational  approaches  and  with  related  mathemati¬ 
cal  and  computational  questions.  It  also  involves  the  questions 
of  redundancy  symmetry  and  the  appropriate  way  tc  describe 
these  fields  for  the  continuous  and  non  continuous  cases, 
which  may  (and  most  often  will!)  occur  for  nonlinear  systems, 
e.g.  flows  with  shocks. 

The  last  question  as  well  as  the  elaboration  on  the  otnei 
points  are  in  the  works  now,  to  be  reported  at  a  later'  date. 

Simple  computational  examples,  where  ' '■•'rnat  i  ve  formu¬ 
lations  and  trial  functions  are  used  for-  the  Laplace  and  i'ri- 
comi  problems  are  reported  in  the  following  chapter. 
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III.  ON  DIFFERENT  MIXED  FINITE 
ELEMENT  APPROXIMATIONS  FOR  SYMMETRIC 
ELLIPTIC  SYSTEMS 


Sara  Yaniv 


Abstract 

Different  finite-dimensional  spaces  are  tried  for  mixed 
finite  element  approximations  for  a  functional  which  has  a 
saddle  point  at  the  solution  of  elliptic  symmetric  linear  system, 
with  a  Lagrange  multiplier.  Calculations  are  carried  out  for  2 
first  order  equations,  (u,v)(x,y)  and  boundary  conditions,  e.g. 
Laplace  and  Tricomi's  problems.  Brezzis'  convergence  condition 
is  found  hard  to  verify  rigorously,  even  for  descretizations  that 
seem  to  work  well.  Preliminary  analysis  is  tried  and  experiments 
conducted  for  bilinear  variations  on  rectangles  for  all  compon¬ 
ents  and  for  bilinear  (u,v)  and  piece-wise  constant  (A)  trial 


functions . 


2 


1.  Introduction  and  variational  formulation 
Consider  the  equation 

Cla)  Ax  +  By  =  f(x,y)  (A,B)  (x,y  ,ij>x  ,<|>  ) 

in  a  domain  il  with  given  boundary  conditions: 

C lb )  j  “  £ 

Assuming 

u  = 

V  =  ^ 

and  A  =  B 

v  u 

then  the  operator  F(x,y,u,v)  =  A  ■+  B  is  a  potential  [1,  p... 

x  y 

35]. 

The  variational  formulation  for  the  problem  is:  find  (u,v;>. ) 
so  that 

(2)  J (u , v; X )  =  I  j[L(x,y,u,v)  +  XCu  -v  ) 3dxdy 
J  Sh  y  x 

-  2  F(x,y)udxdy 
n- 

x 

F(x,y)  =  f(£,y)d£ 

is  stationary,  for  all  functions  (u,v)  £  VxV  and  x  £  W ,  where 
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e'-ery  (u,v)  €  VxV  satisfying  boundary  conditions: 


%  'X. 

udx  +  vdy  =  0 . 


Lvx,y,u,v)  is  the  Lagrangian  of  the  problem,  i.e. 


L  =  A 
u 

L  =  B 
v 

and  x(x,y)  is  a  Lagrange  multiplier  which  is  the  stream  fun 
of  the  problem  [1,  page  3S]. 

The  variation  of  J(u,v;X)  for  fixed  X  gives  the  rci 
ing  weak  formulation: 

f 

(  ■O  5J(u,v;X)  =  [L  6u  +  L  6v  +  X(6u.,-6v  )  ]dxdy  =  0 

•  u  y  x 

Adding  the  equation 

u  -  v  =  0 

y  x 

in.  weak  formulation: 

r  r 

|  q(Uy~v^)dxdy  =  0  for  Vq£W , 

we  get  a  saddle-po.i.nt  problem: 

f  r  9U-.Bv, 

(4a)  I  [L  u  +  Ly  v  +  *(  —  -  — )]dxdy  =  0 
1  ft' 

v  (uL  ,v1>  e  v  x  v 
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and 

(4b)  |  qC|^  -  j^Odxdy  =  0  Vq£W. 

The  functions  u,v  and  A  which  satisfy  (la)  and  (lb)  are 
the  solution  of  C4a)  and  CHb). 

For  linear  problems: 


5 


This  formulation  gives  a  saddle-point  problem  [2],  which  is 
determined  as  follows: 

Given  f  €  V'  ,  g  £  W'  ,  find  u  €  V ,  \p  €.  W  such  that: 

a(u,v)  +  bCv,i/0  =  f  (v)  Vv€V 

(5) 

b  ( u  ,<t>  )  =  g(  <p)  V<t>£W 

Brezzi's  theorem  states  the  following  conditions  for 
existence  and  uniqueness  of  the  solution: 

Suppose  a  and  b  are  bounded  and 

C6a)  inf  s up  |a(u,v)|  5  y  >  o 

u€z  v£z 
ilu  |!  v  =  1  ||v  11^  =  1 

z  :  {u£V/b(u,<t>)  =  0  V<J> GW } 

t6b)  sup  * }  I  5  k  114;  U w  V*£W,  k  >  0 

v€  V  IIVIIV  W 

then  the  solution  of  (5)  is  unique 
and 

«u#v  +  ^  «w  s  C(  If  lv  ,  +  Jg#w>)  . 

It  is  easy  to  verify  that  problem  C:a),  (4b)  satisfies 
(6a),  (6b),  (6c),  hence  has  a  unique  solution  which  is  the  only 
solution  of  (la),  (lb). 


2.  Approximation 
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In  order  to  approximate  the  solution  of  problem  (4a),  (4b) 
using  the  saddle-point  weak  formulation  or  the  saddle-point 
variational  principle,  we  use  finite-dimensional  spaces  c  V, 

c  W  satisfying: 

(7a)  inf  sup  |a(u,v)j  >,  t  z  0  ,  x  independent  of  h 
u£Zh  v6Zh 
nui  v=  i  «v»v=i 

(7b)  sup  |  a(  u,  v)  I  >0  V  0  *  u  E  Z, 

vEZ,  h 

h 

Z^  =  (u£Vh;  b(u,4>)  =  0  V(j!ewh) 

<7c>  vev  'Lb.".tv)|  J  ‘  '*'h 
h 

v^ewh 

£  >  0 ,  independent  of  h. 

The  approximated  problem  is: 

a(uh,v)  +  b(v,<|>h)  =  f  ( v )  VveVh 

(8) 

b ( uh , ^ )  =  g(4>)  V4>€Wh  . 

Under  hypothesis  (7a),  (7b),  (7c)  problem  (8)  has  a  unique 
solution  [2]  and: 

|u  -  Vy  +B4-  -  Vw*  c  inf  (|U-X,V+|*-6|W 

xtvh 

66Wh 
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3.  Examples 

We  used  the  variational  principle  (2)  to  solve,  approximate! 
the  Laplace  equation  and  the  Tricomi  equation  in  an  elliptic  doma 

Ci)  The  Dirichlet  problem  for  the  Laplace  equation 


>  +  $  =0 
xx  yy 


3ft 


=  f(x,y) 


The  variational  functional  is: 


C9)  J(u , v ; X )  = 


[u2  +  v2  +  x(u  -v  )]dxdy 
y  x 


for  (u,v)  €  VxV  ,  X  €  W 


VxV  ^  { C  u , v ) / 


2  2  2 

(u  +  v  +(u  -v  )  ]dxdy  <  °° ,  udx  +  vdy 
y  x 


=  di 


3  £3 


w  =  {q  | 


q2dxdy  <  «  ,  q(xQ ,yQ)  =  0} 


for  this  problem  conditions  (6a),  (6b),  (6c)  are  fulfilled, 
hence  the  problem  has  a  unique  solution. 

Let  ft  be  the  rectangle: 


ft  =  { (x,y )/  -1  *  x  s  1,  -1  s  y  s  0} . 
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We  have  divided  fl  into  rectangles  of  size  h  and  chosen  dif¬ 
ferent  finite  dimensional  spaces  for  the  trial  functions 

a)  The  first  attempt  was  to  use  the  bilinear  trial  func¬ 
tions  for  (u,v)  and  X,  for  these  spaces  Brezzi’s  conditions  are 
not,  necessarily  satisfied;  this  is  concluded  from  the  instability 
of  the  numerical  solution,  in  some  of  the  cases  tried, 

b)  We  tried  to  take  other  finite  dimensional  spaces  so  that 

(7a)  ,  (?b),  (7c)  will  be  satisfied.  Using  the  same  finite-element 
discretization  of  as  in  a) ,  approximating  u  and  v  by  bilin¬ 

ear  trial  functions  and  piece-wise  constant  functions  for  X 
(intuitively,  this  may  help  to  get  rid  of  the  4  constants  and 
leave  us  with  the  only  arbitrary  constant  of  the  problem  which  is 
qCx0,y0)  =  0). 

For  this  approximation  the  solution  converged  to  the  analytic 
s Glut  ion ; 

X  .  X  X  «•  l 

u  -  e  siny,  v  -  e  cosy,  A  =  -2e  cosy  +  2/e  as  shown  in  the 


following  table: 
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L^-errcr 


h  =  0.25 


L--  error 


h  =  0.125 


L~-erroi 


h  =  0.0625 


numerical  rate 
of  convergence 


0.0.415 


0. 016  3 


0.053 


0.0194 


0.184 


0.097 


1.  35 


1.45 


0.92 


(The  rates  of  convergence  for  all  of  the  problems  we  solved 
were  the  same). 

In  [3,  page  75-77]  the  authors  show  that  for  the  Stokes 
equation  (for  a  mixed  variational  formulation)  quadratic  trial 
functions  for  u  and  v  and  piece-wise  constant  function  for  k 
for  triangular  elements  are  permitted,  and  then  conditions  (7a), 

( V b ’ ,  (7c)  are  satisfied. 

r)  'tie.  t  f-i.ed  to  solve  the  variational  principle:  find  (u,v;') 
so  that 


10 


J(u,v; A) 


A  u  +  A  v]dxdy  + 
y  x  “ 


+ 


3fl  A(fxdX  +  fydy) 


is  stationary,  for: 


(u,v)  £  VxV  ,  A  e  w 

where 

VxV  =  {(u,v)/[  f(u2  +  v2)dxdy  <  °°} 

w  =  {q/j  l(qy  +  qx)dxdy  <  "  ’  tJ(xo>yo)  =  0}> 

as  an  equivalent  problem  for  the  Laplace  equation. 

For  this  f  oirniu'lation  we  used  the  same  element  discretiza 
tion,  but  took  piece-wise  constant  trial  functions  for  u  and 
and  bilinear  trial  functions  for  A.  The  numerical  solutions 
unstable,  we  received  2  independent  solutions  for  A  and  the 
solution  for  u  and  v  did  not  converge. 


v 

was 
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ii)  The  Dirichlet  problem  for  the  Tricomi  equation  in  an  ellip¬ 


tic  domain 


yi  -  it  =  0 
J  XX  *yy 


=  f(x,y) 


Q  =  {(x,y)/  - 1  *  x  $  1 ,  -1  $  y  §  0 } 


The  variational  functional  is: 


f  f  o  o 

CIO)  J(u,v;A)  =  Cyu^  -  v  +  A(u  -  v  )]dxdy 

Jo  j  y  x 


for  (u,v)  €  VxV  ,  A  €  W 


VxV  =  {(u,v)/  [yu2  -  v2  +  (u  -  v  )2]dxdy  <  «® , 
n.J  y  x 


udx  +  vdy  =■  df) 
30 


2 

W  =  {q |  (q  dxdy  <  °°  ,  q(x0>y0)  =  0} 


conditions  (6a),  (6b),  (6c)  are  satisfied,  hence  there  exist  a 
unique  solution. 

We  followed  the  numerical  procedures  a)  and  b)  used  for'  the 


Laplace  equation.  The  approximated  solutions  behave  the  same 
Cl,  page  73]. 
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For  procedure  a)  the  solution  for  A  was  unstable  dependin 
on  4  arbitrary  additive  constants  while  for  u  and  v  the 
approximation  was  good. 

For  Procedure  b)  the  approximated  solution  was  stable  and 
converging  for  u,  v  and  A  as  shown  in  the  following  table: 


The  analytic  solution  in  this  case  is: 


u 


v 


sinhX<yT  fon»l)3n'<3n  -? )<3n-3) .  ..4.3 
00  v3n 

coshx( l  +  £  3n( 3n-2 7T 3n-3 ) . . .4  71  > 


»  3n 

A  =  2sxnhx(l+£  3n  (  3ri-  2  M  3  n-  3  5  .  '.'.'n  .3  > 


The  rate  of  convergence  is  the  same  for  all  the  case 
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4.  Remarks 

Equation  (la)  is  equivalent  to  the  first  order  system: 


0 
0 

0 

for  the  Laplace  equation 

for  the  Tricomi  equation 

It  is  obvious  that  the  solution  for  X  has  more  derivative 
than  u  and  v. 

The  trial  functions  used  in  procedures  a)  and  b)  do  not 
satisfy  this  feature  while  those  of  c)  satisfy  it.  But  the  only 
procedure  which  gave  a  converging  solution  is  b). 

The  only  conditions  which  insure  convergence  of  the  approx! 
mated  solutions  are  those  of  Brezzi  (7a),  (7b),  (7c),  (which  are 
quite  difficult  to  show  in  the  finite-dimensional  problem) . 


L  -  X 
u  y 


(11) 


L  +  X 

V  X 


where 


and 


u  -  v 

y  x 


L  =  2u 
u 


L  =  2v 
v 


Lu  =  2yu 


Lv  =  -2v 
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Summary  and  Concluding  Remarks 

A  summary  of  results  u  finite  elements  uul 
non-uni  formJy  elliptic  pa:  is  are  riven  in  !.  j  ’ .  A 
gion  is  descretized  into  rectangles  and  a  e :  line 
assumed  for  both  field  variables  lu,v)  ana  the  h 
X  appearing  in  the  functional: 

C12)  J(u,v;X)  *-  j  J  [L  +  >.  tu  -  v  )  1  dxdv 

J y  x 

to  be  made  stationary  at  the  solution  l unctions 
The  values  of  X  at  one  c.anei  <  i  -  J,  ; 
mine  its  values  at  even  points  only  ana  tne  set 
the  odd  points  consist  of  a  separate  oyster:  ,  It  t 
cep  table  behavior  of  the  overall  solution  (a  phe 
other  computational  contexts).  To  coup  re  the  equ 
points  the  values  of  \  at  the  4  points  of  the 
had  to  be  predetermined  (e.g.  by  a  Taylor's  ex pa 
CQ,C)  point  and  the  relation  to  the  values  1  u, 
rectangle).  The  procedure  is  not  considered  sati 
it  yields  acceptable  results.  On  the  hyperlc-lit 
in te rrr.it tency  causing  unstability  has  been  cbs»-> 
Following  a  series  of  lectures  by  J.  Osborn 
variation  was  retained  for-  (u,v)  but  the  tent  q. 
replaced  by  piece-wise  constants  on  rectangles, 
occurs  and  the  calculations  are  completely  stall 
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The  price  is  paid  in  a  lower  accuracy,  but  the  scheme  is  never¬ 
theless  considered  feasible. 

The  functional: 


J( u, v; X ) 


[L  +  uX  -  vX  ]dxdy  +  X  (  f  dx  +  f  ,cy) 
y  x  j  x  y 

aa 


obtained  from  (12)  by  integrating  the  2—  term  by  parts,  is 
also  stationary  at  the  solution  to  the  Laplace  system.  The  same 
des critization  and  the  corresponding  trial  spaces  yield  unstable 
scheme.  This  is  somewhat  unexpected,  because  the  relative  differ¬ 
entiability  required  for  (u,v)  and  X  in  (11)  fits  better  the 
analytic  relation  than  the  one  in  b).  The  stability  analysis 
for  this  case  (via  Brezzi's  coercivity  condition  for  the  approxi¬ 
mation  spaces)  is  difficult,  and  we  have  not  been  able  t.  c^.  . 
it  to  a  successful  conclusion. 


| 


! 
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IV.  FIN  IT  I  IIEP-  1LI.CI.  APPROXIMATIONS  FOR  THL 

SOLUTION  OF  TH-.  TPJCOMI  EQUATION  IN  A  MIXED 
ELLIPTIC-HYPERBOLIC  REGION 
by 

David  Levin 


Abstract 

Difference  schemes  for  the  solution  of  the  Tricon i  pi 
are  derived,  in  both  the  elliptic  and  hyperbolic  domains  : 
schemes  are  constructed  so  to  be  exact  for  several  pcivnoi.ds. 
solutions  of  the  Tricomi  equation.  The  method  presented  i;: 
adaptable  to  non-luv-M  ornr  i  ns.  ar.d  tc  non-standard  meshes. 
A  high  accuracy  is  do  men  at  rj'  eu  lor  .  everal  Tricon,  i  boundary 


conditions . 


1. 


Introduction 


The  Tricomi  equation 

>';xx  -  Jyy  =  C  <  — > 

is  elliptic  for  y  <  Cl  and  hyperbolic  for  y  >  0  with  char 
acteristics  defined  by: 

f«  =  x 

l 

n  r  x 

Tricomi  [6]  showed  that  (1.1)  has  a  unique  solution  in  a 
domain  D  =  I)+  Cl  D  bounded  by  a  simple  arc  F  in  the  ellip 
tic  domain  with  endpoints  on  the  x-axis  at  x  =  A  and  x  =  £ 
and  by  the  two  characteristics  r,  and  r0  through  these 
endpoints  (figure  1.1). 


3/2 


3/: 


(1.2) 


Figure 


Problem  (1.1)  is  known  to  be  well  posed  in  D  with  Dirich 

let  boundary  conditions,  i.e.  ,  ■;  given  on  r~  anc  cn 

Fq  h  { y  =  0 ,  A  f  x  E }  .  In  I  problem  (1.1)  is  well  posed 

with  either  Cauchy  conditions  -  c  and  <pv  given  on  or 

Goursat  conditions  -  $  giver,  on  r.  and  on  either  r,  or  r„ 

J  1  l 

For  solving  the  Trieonsi  problem  in  Z  one  should  match  scheir.es 
in  D+  with  schemes  in  D  ,  e.g.  .  Filipov  [1]  usee  a  i'ive-pcint 
formula  in  D  and  a  four-point  roimula  in  E+  matched  by  a 
difference  equation  corresponding  to  *  =  0  along  the 

"parabolic  line"  ,  on  a  mesh  as  in  figure  2.1 

Figure  1.2 


Filipov  also  proved  that  his  scheme  is  regular  and  converging 
The  change  in  the  structure  of  the  mesh  across  the  parabolic 
line  makes  it  difficult  to  construct  higher  order  schemes  cl 


this  kind  in  the  usual  manner. 
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Another  strategy  for  solving  the  Tricomi  problem  is  that 
of  Vincenti  and  Wagoner  [7]  who  reduced  the  problem  to  a  pure 
elliptic  problem.  However,  the  boundary  conditions  on  r  , 
which  are  obtained  by  projecting  the  given  conditions  on  r ^  > 
are  quite  complicated  and  cause  great  numerical  difficulties. 
Several  authors,  [2],  [5]  have  used  expansions  in  terms  of 
certain  particular  solutions  of  (1.1),  and  this  method  proved 
to  be  quite  effective  for  cases  of  very  smooth  boundary  condi 
tions . 

The  method  presented  in  this  work  combines  somehow  the 
motifs  of  the  above  three  strategies;  a  local  expansion  in 
terms  of  particular  solutions  of  the  Tricomi  equation  are  use 
to  produce  high  order  difference  schemes  for  the  Tricomi 
problem.  Using  these  schemes  the  problem  is  reduced  to  an 
elliptic  problem  in  D  with  certain  boundary  conditions  on 
y  =  0. 

The  method  used  here  for  producing  the  difference 
scheme  has  recently  been  found  to  be  useful  for  solving  both 
Cauchy  and  Goursat  problems  in  D+ .  It  is  the  versatility  of 
this  method  which  enables  us  to  obtain  specially  structured 
difference  schemes  of  high  order  for  the  desired  matching 
along  the  parabolic  line. 

In  section  4  we  describe  some  numerical  experiments  with 

3 

the  suggested  procedure,  exhibiting  a  global  9(h  )  accuracy. 
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2.  The  expansion  method  for  producing  difference  schemes 
Let  z^  =  (x^,y^)  i  =  +1  be  N  +  1  adja¬ 

cent  mesh  points  in  a  grid  covering  the  domain  D.  We  look 
for  a  difference  approximation  for  the  Tricomi  equation 

(1.1)  based  upon  {z.}^+^.  The  usual  way  of  obtaining  dif- 

1  i  =  1 

ference  schemes  is  by  expanding  the  approximation  u(x,y) 
in  power  series,  around  z^+-^  for  instance,  then  form  a 
linear  combination  of  u(  )  , .  .  .  ,u(  z^+  ^  )  ,  using  (1.1),  to 
get  the  desired  degree  of  accuracy.  This  process  can  be 
simplified  by  using  an  expansion  in  terms  of  the  polynomial 
solutions  of  (1.1).  These  can  be  found  by  expressing 
4>(x,y)  as  a  double  power  series  expansion  around  (0,0)  and 
comparing  the  expansions  of  y<f>xx  with  that  of  <f>yy  .  The 
result  is  that  all  the  polynomial  solutions  of  the  Tricomi 
equation  can  be  written  as 


P2M+l(x’y) 

■P2M+2^x,y^ 


.  . 

r  M-2 i  3i 
Z  c.x  y 

i  =  0 


[M] 

?  .  M-2 i  3i+ 1 

I  d.x  y 

i=0 


M  =  0,1,2,.  .  . 


(2.1) 


where  Cq  =  d^  = 


1  and 
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r „  .  (M-2i) (M-2i-l) 

i+1  '  (3i+3) C3i+2)  ci 


^  _  (M-2i) (M-2i-l)  , 

i  +  1  "  (3i  +  4) ( 3i+  3 )  ai 


i  =  1,..  . ,lyj. 


(2 .2) 


LEMMA. 

Let  u(x,y)  be  an  analytic  solution  of  (1.1)  in  a  neigh¬ 
bourhood  S  of  (xQ,y0)  then 


2M+1  M 

u(x,y)  =  £  a.P.(x,y)  +  o((max(h,k)  )  V  (x,y)  £  S  (2.3) 

3=1  3  3 


where  h  =  |x-xnj  and  k  =  |y-yn|.  Also,  for  yn  =  0  and 


k  *  h' 


2K+1  M 

u(x,y)  =  l  b.P.(x,y)  +  °(hn) 
j  =  l  3  3 


(2.4) 


2M  +  2 

l  b . P . ( x  ,y  )  +  °(h  ). 

3=1  3  3 


The  proof  is  straightforward. 

Using  the  above  lemma  it  is  clear  that  if  we  find  a 

N 

scheme  which  is  accurate  for  (Pj)j_^,  i.e. 


I  aipj(xi,yi)  =  Pj(xN+1>yN+1)  >  3 

1=  1  J  J 


=  1,2,. . .  ,N,  (2.5) 


then  this  scheme  has  a  truncation  error  of  order  [  y]  There- 
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fore  we  simply  use  (2.5)  as  the  defining  equations  for  our 
schemes.  This  method  of  obtaining  difference  schemes  is 
convenient  to  use  for  any  given  distribution  of  mesh  points, 
provided  that  the  system  (2.5)  has  a  solution.  The  method  has 
recently  been  used  successfully  for  solving  the  Cauchy  and 
Goursat  problems  in  D+  [3],  However,  the  schemes  presented 
in  [3]  are  not  in  a  suitable  form  to  be  used  for  solving  the 
mixed  problem  in  D  since  they  cannot  be  matched  nicely  with 
schemes  in  the  elliptic  domain  D  .  In  order  to  obtain  suitable 
schemes  we  consider  the  "discrete  Cauchy  problem"  in  D+ ,  i.e. 
solving  (1.1)  with  the  two-level  conditions 

/  4>  ( x  , 0  )  =  f-^(x) 

j  AfxjB  (2.6) 

U  (x  ,-<5  )  =  f  2  (x ) 

.  •  +  . 

This  problem  is  not  well  posed  in  D  m  the  usual  sense. 

However,  if  we  recall  that  <J>  should  be  a  solution  of  the 
mixed  problem  in  D,  ♦  should  be  in  in  a  neighbourhood  of 
y  =  0  for  any  A  <  x  <  B.  Therefore,  $  (x,0  )  =  $y(x,0  )  and 

+ 

as  6  -*■  0  the  problem  (2.6)  turns  to  be  a  Cauchy  problem  in  D  . 
For  a  fixed  6  we  should  also  give  the  additional  boundary 
conditions : 

l'*(A,y)  =  g1(y) 

<  -6  <  y  <  0  (2.?) 

U(B,y)  -  g?(y) 
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to  make  the  problem  well  posed  in  A  s  x  *  B  -6  %  y  t  C ,  and 
hence  also  in  D+ . 

For  a  small  6  it  is  expected  that  the  domain  of  influ¬ 
ence  of  the  boundary  conditions  (2.6)  is  similar  to  that  of 
the  Cauchy  conditions.  Therefore,  we  consider  a  mesh  defined 
by  characteristic  lines  in  D+  as  in  Tl]  and  [3],  lienee,  wc- 
look  for  a  scheme  of  the  form 

M 

l 1Cai*(Xi»ym)  +  aM+i*Cxi’ym-l)]  =  ,ym+l5  (2‘8) 

i=  1 


m  -  0  ,1 ,  .  .  .  ,2n-1  ,  where 


y 


-l 


-  6 


6  t  (|h) 


2 

3 


yc  =  o, 

i  I 

vm<1  =  4h  *  y2J 3  "  =  0.1 . 2S-! 

and 

xi  =  x0  +  h(fcl  "  i  =  1.2,.  ..,M 


whe  re  h 


B-A 

7K 


as  in  f igure  2.1. 


with  N  =  2M,  i.e.  by  the  system  ol  2M  linear  eqi 


J  [  a  .  P  .  ( x  -  ,y  )  +  a .  (*:•: .  ,  v  ,  )  !  ~  P  .  C  xn  , 
l  :  i  rr.  i  ;  :  r.-i  3  O’ 


Given  the  "discrete  C<v;chv  c.r.^ 
scheme  (2.6)  for  a  =  1  (  aC.  cumin -  hut  (2.9)  has 
one  can  get  an  upproxin.ati  •  r,  : ■  tn<-  level  v  = 
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to  move  upward  to  at  the  top  of  D  .  Before  proceeding 

to  the  mixed  problem  we  present  seme  numerical  results  with 
the  above  scheme  for  the  "discrete  Cauchy  problem";  we 
examined  the  case  M  =  3,  i.e,  the  scheme  is  built  to  be  exact 
for  P,,...,Pg.  However,  it  so  turns  out  that  the  system  (2.9) 
is  singular.  Therefore,  we  could  make  the  scheme  accurate 
for  P and  Pg  as  well,  and  thus,  by  (2.3)  and  (2.4)  the  re- 

suiting  scheme  is  of  order  o(h°  °)  near  y  =  0  and  of  order 

3  3+2/3 

°(y  ,  ,-y  )  )  for  mil.  A  global  0(h'-  )  accuracy  has 

been  detected  in  a  series  of  numerical  experiments.  An  example 

of  this  is  shown  in  the  following  table  where  an  analytic 

solution  cf  the  Tricomi  equation, 

“  3n+l 

t(x,y)  =  cosh  x(y  -  [  Hn"i )  3'n  •  (In- 2  TTTu  "3  }  (2  —  } 

n-  1 

is  taken  as  a  test  function.  The  table  shows  the  maximum  ab¬ 
solute  value  of  the  error  obtained  by  using  the  suggested 
scheme  with  h  =  i,  where  the  "discrete  Cauchy  condi¬ 

tions  (2.6)  are  given  with  6  =  j. 

Table  2.1 


1 

1 

T_ 

1 

h 

2 

4 

*8 

16 

max | error | 

. 67E-4 

. 57E-5 

.  39E-6 

25E-7 

1  a  - 


3.  The  scheme  for  the  mixed  problem 

We  now  consider  the  Tricomi  problem  in  the  mixed  domain 
D  =  D+  U  D~  with  boundary  conditions  on  r  and  on  r , . 
Assuming  that  r”  is  rectangular  we  can  use  a  rectangular 
mesh  in  D-  and  5  or  9  point  schemes  around  any  internal  mesh 
point  in  D  .  The  scheme  (2.8)  for  m  =  0  can  be  regarded  as 
difference  equation  for  the  mesh  points  on  y  =  0.  However’, 
when  we  move  on  upwards  we  find  that  we  still  miss  some  diff¬ 
erence  equations  to  stand  for  the  unknowns  at  the  mesh  points 
on  the  free  boundary  T^.  Therefore,  some  additional  schemes 
should  be  introduced. 

Let  us  denote  by  cp^  the  vector  of  the  values  of  the  ap¬ 
proximation  at  the  mesh  points  on  y  =  0  and  by  ip_^  the  vec- 

^  ,  .  .  ,  , 2N( M-l ) - 1 

tor  of  the  values  on  y  =  -o,  a.e.  <pQ  =  an- 

».  =  where  1=1 

1  11  i=l 


r«Poi  =  $<A  +  ippp  0) 


l'Pli  =  <>(A  +  _(5) 


i  =  1  , .  .  .  , 2N ( M- 1 ) - 1 . 


(  2  .  J  ) 


Successive  use  of  the  schemes  (2.8)  for  m  =  C,l,...,2N-2 
finally  yields  a  linear  relation  between  the  approximation  at 
the  mesh  points  on  and  iPg  and  tp  ^  in  the  form: 


*(A h,y  )  =  Si  m  +  Tl  ip 


-  qT 


T 

•r  J 


m 


=  ,2U 


t 
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where  S  and  T  are  ( 2N ( M- 1 ) -1 ) -dimensional  vectors. 

Since  4>  is  given  on  F^  we  can  consider  (3.2)  as  2N  equa¬ 
tions  relating  the  unknown  coefficients  in  tp^  and  (p  . 

Yet  we  are  short  of  2N(M-2)-l  equations  to  complete  the  sys¬ 
tem.  To  get  these  equations  we  use  some  additional  schemes 
for  solving  the  "discrete  Cauchy  problem"  which  result  in 
relations  of  the  form  (3.2)  between  tpn  ,  tp  ^  and  bcunaary 
values  at  intermediate  points  on  r^.  This  can  be  done  for 
any  M  and  with  no  less  accuracy  than  in  (3.2).  To  demonstrate 
it  we  present  the  treatment  for  the  case  M  =  3; 

For  M  =  3  we  are  short  of  2N-1  equations  and  we  could 
have  just  the  right  number  of  equations  if  we  could  get  re¬ 
lations  of  the  form  (3.2)  for  the  boundary  values  <t>  ( ,y  , 
where 


m+H 


m+i 


A  +  2iih 


3  2 

,  3,  .  2.3 

<8*  *  V 


m  =  1,  .  . .  , 2N-1 


(3.3) 


Such  relations  can  be  obtained  by  using  additional  4th  order 
schemes  of  the  form 


4 

.Z  ,y  }  +  “4  +  i^^.y  3)]  =  *(x  3  ,y  3>-  <3.4) 

i=l  x  m+T  m+  — 

/.  c 


Again  the  coefficients  of  these  schemes  are  obtained  by  sys- 
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terns  of  equations  of  the  form  (2.5)  and  here  we  could  make 
the  schemes  be  accurate  for  P^,...,Pg. 

In  the  end  we  have  the  2N(M-1)-1  relations 

<KA  +  yn>  =  S^Q  +  T^_x  n  =  2  .  .  ,4N,  (3.5) 

2  2  2 

where  S  and  T  are  given  ( 2N(M-l)-l)-dimensional  vectors, 
n  n  6 

2  2 

These  4N-1  relations  (for  M  =  3)  are  combined  with  the  ordinary 
scheme  for  the  internal  points  in  D  to  give  a  full  system  of 
equations  for  all  the  unknowns  in  D  and  on  the  parabolic  line 
y  =  0.  After  solving  this  system  we  can  use  the  vectors  and 

n  :  2  , .  . .  ,  4N  to  produce  an  approximation  at  all  the  mesh  point 

T 

.  + 
in  D  . 
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4.  Numerical  experiments 

In  this  section  we  present  some  numerical  results  of 
applying  the  new  schemes  for  several  Tricomi  problems.  We 
used  the  presented  schemes  (with  M  =  3)  in  the  hyperbolic 
domain  combined  with  a  simple  5  and  9-point  formula  in  the 
elliptic  domain.  The  3-point  formula  is  also  obtained  by  a 
procedure  of  tire  type  (2.5),  i.e.,  by  dem.anding  that  the 

scheme  is  accurate  for  polynomial  solutions  of  the  Tricomi 

4  - 

equation.  In  that  way  we  obtain  a  local  °(6  )  scheme  in  b 

4 

which,  experimentally,  proved  to  be  a  global  0(6  )  where  a 

square  mesh  of  size  6  is  usee  in  D- .  In  order  to  match  the 

hyperbolic  and  elliptic  schemes  we  chose  6  =  y. 

We  consider  boundary  conditions  given  on  the  lines 

r  _  =  { x  =  1 1 ,  - 1  y  £  0 }  and  { y  =  - 1 ,  - 1  s  x  s  1 }  ir.  *. :  • 

elliptic  domain  and  on  the  characteristic  line 

2  2 

2  2  3  r 

r i  -  { x  -  -jy  =  -1,  0  <  y  5'  (^O'}  in  the  hyperbolic  domain . 
Such  conditions  define  a  unique  solution  in  the  domain  bounded 
by  the  above  lines  and  by  the  characteristic 

3_  ?_ 

-  (x  +  -|y2  =  1,  0  <  y  <  (y)3}. 

The  computational  aspects  of  the  method  and  some  addition  a  1 
numerical  results  are  to  be  described  in  a  separate  rep  or4  !  u  .! 
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Table  3.1 

Maximum  error  in  D  =  D+  U  D  for  the  Tricomi  problem 
with  the  values  of  the  analytic  solution  (2.10)  of  the  Tri¬ 
comi  equation  as  boundary  values  on  U  r  .  A  scheme  with 

M  =  3  is  used  in  the  hyperbolic  domain  and  5  and  9-point 
formula  are  used  in  the  elliptic  domain. 


hl'l  a 

A>:  =  Ay  =  j  If  g  r6 


5-point  formula 

0. 47E-2 

0 . 14E-2 

0  .  37E-3 

9-point  formula 

0 . 44E-3 

0.5SE-4 

0  .  “81-5 
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V.  DIFFERENCE  SCHEMES  FOR  THE  SOLUTION 
PE  TRICOMJ '  S  ECU  AT  IP!-;  i.N  A  MIXED  HECIO'.' 


Frieda  Loinger 

Abstract 

The  Tricomi  problem  in  a  mixed  region  is  solved  via  a 
numerical  projection  of  the  hyperbolic  boundary  conditions 
onto  the  parabolic  line,  thus  coupling  the  two  regions.  Differ¬ 
ence  schemes,  exact  for  polynomials  up  to  a  certain  degree, 
are  used;  high  accuracy  is  achieved  in  all  the  examples  compu- 


t 
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1.  Introduction 


Tricomi's  equation: 


ytp 

;;x 


ip  =0 

yy 


is  elliptic  for  y  <  0,  parabolic  for  y  =  0  and  hyperbolic  for  y  >  0 
We  look  for  a  solution  in  a  mixed  domain  D  =  D  U  D+  where  D_  is 
an  elliptic  rectangle  bounded  by  x  =  -1,  y  =  -1,  x  =  +1,  y  =  0  and 
D+  the  hyperbolic  curved  triangle  bounded  by  y  =  0  and  the  two 
characteristics  1’^  and  r^: 


1  • 


-1  =  n  =  X  -  rry 


2  3/2 


-1  £  x  £  0 


,  .  .  2  3/2 

1  =  C  =  x  +  3  y 


0  £  X  £  1 


The  problem  is  assume  well-posed  with  appropriate  boundary  con¬ 
ditions  given  on  x  =  -1 ,  y  =  -1 ,  x  -  +1  and  one  of  the  character¬ 
istics,  say  r,  • 
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2 .  Notation 

For  the  numerical  solution  of  the  problem  we  divide  the  region 
into  the  following  elements  (see  figure  1):  rectangles  in  the  ellip¬ 
tic  domain  along  the  lines  x  =  const,  y  =  const  and  isoparametric 
triangular  elements  in  the  hyperbolic  domain,  with  the  following 
enumeration : 

j  -  the  no.  of  the  row 
j  =  -M,  -M  +  1  ,..._, 0,1, ...  ,N 
j  <  0  elliptic  region 
j  =  0  parabolic  line 
j  >  0  hyperbolic  region 
*  is  a  mesh  point 

no.  of  points  on  the  parabolic  line  :  2N  +  1 


n .  - 

3 

the  no. 

.  of 

point 

s  in 

the  j-th-row 

xi  +  l 

-  x .  = 

l 

Ax 

=  1/N 

Vi 

in  the  grid 

yj+l 

-  v3  = 

Ay 

=  1/M 

for 

j  <  0 

y  j  +  1 

-  y3  = 

Ay. 

=  (1. 

5  Ax  ( ; 

.±1,,2/3  ,,  ,  A  ..2/3  f 

)  +1) )  -  (1.5  Ax *3)  for 

j  i  0 


mesh  points  in  the  elliptic  part  including  the  parabolic  line: 


3 


(xi ,yj )  i  =  -N,  -N  +  1 ,  .  .  .  , 

x .  =  iAx  v • 

i  y  3 

mesh  point  in  the  hyperbolic  part: 

(x.  • ,y . )  j  =  1,2,. .  .  ,N 
-1  5  J  J 

y.  -  ( 1 . 5 Ax j ) - ^ ^  ;  x 
and  the  numerical  solution  at  the 
)  =  ip(x  .  ,y. )  i  = 

X  >  J  J 


1  >1. • • • ,N  ;  j 

=  jAy 

i  =  -N, . . 

•  =  ( i+ j ) Ax 

>  1 

-th-row : 

N ,  -N  +  1 , .  . 
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3 .  Projection  of  the  boundary  condition  along  r-,  on  the  para¬ 
bolic  line 

The  hyperbolic  problem  with  either  Goursat  conditions  (<p  given 
on  ri  and  y  =  0 )  or  Cauchy  conditions  C  tp  and  tp  given  on  y  =  0)  is 
treated  in  [1]  (D.  Levin)  using  high  accuracy  difference  schemes. 
The  difference  scheme  for  (x^,  y^  +  ^)  is  based  on  the  values  of 

at  the  points  ( x  ^  +  Ax,  y..  )  ,  (x^jy^)  when  Cauchy  conditions 
are  given  (figure  2). 

The  coefficients  are  determined  by  the  demand  that  the  scheme 
is  accurate  for  6  polynomial  solutions  of  the  Tricomi  equation: 
l,x,y,xy,3x2  +  y3,x3  +  xy3  .  Numerical  experiments  showed  that 
the  schemes  obtained  are  accurate  for  S  basic  functions: 
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l,x,y,xy,3x  +  y  ,x'  +  xy  ,6x  y  +  y  ,2x  y  +  xy  . 


We  now  look  for  a  solution  in  the  mixed  region. 

An  analytic  connection  between  the  boundary  conditions  on  r  ^  and 
the  oarabolic  line  is  known: 

Bitsadze  [2J: 


t  ( x )  =  ip(  x  ,  0  )  = 


2  n  y. 


,5/6 

dx 


*  ijj  ( t  /  2  ) 

0  t2/3(x-t)1/r 


dt 


+  Y 


» (  t ) 


0  (x-t) 


1/3 


dt 


5 


where  v(x)  =  ip  (x,0) 

y 

ip  (  x  ,  0  )  =  <p(x,y(x))  on  T, 

The  connection  on  the  parabolic  line  enables  one  to  find  the 
solution  of  Tricomi's  equation  in  the  elliptic  region  and  then 
in  the  hyperbolic  problem  can  be  solved  with  either  as  a  Goursat 
or  a  Cauchy  boundary  conditions. 

This  report  presents  a  numerical  process  for  finding  a  con¬ 
nection  between  the  hyperbolic  and  elliptic  region,  without  using 
<p^  (see  also  [3]).  Alternatively,  we  look  for  a  connection 
between  the  values  on  and  <p  ^  ^ ,  (p^  ^  . 

Assuming  the  elliptic  problem  already  solved,  i.e. 

<p^-^  =  (vp^  \  +  1  >  *  3  =  known,  where  is  the 

solution  on  the  parabolic  line.  We  proceed  from  row  (j)  to 
(j  +1)  (j  >  1)  by  choosing  6  basic  functions  as  before  (in  the 
Cauchy  problem),  but  instead  of  using  values  of  cp,iPy  on  y  =  y.  as 
it  is  done  in  [1]  for  the  Cauchy  problem,  we  take  an  additional  row 

y.  and  the  difference  scheme  will  be  based  on  the  values  of 

J  J-l 

cp  at  the  points  (figure  3):  (x^,y_.+^)  ;  (x.  +  Ax  ,y^  )  ,  (x^  ,y .  )  ; 

(x.  i  Ax,y^_1),(xi,y^_1)  ; 

and  we  require  that  the  difference  scheme  is  exact  for  the  6  basic 
functions.  numerical  experiments  show  that'  these  7-part  differ- 


o 


ence  scheme  using  6  basic  functions  are  exact  for  8  polynomials. 
Note :  if  j  =  0 ,  y j  =  0,  '  ~  Ay ,  i.e.  the  first  row  of  the 

elliptic  region  is  taken  . 

The  difference  scheme  can  be  expressed  by 

3 

+  s*  ,p(x4»yj-i)  =  ,p(xi»yj+i) 

x^  =  x^  +  ( £- 1 ) Ax  i  =  -N,...,n. 


The  coefficients  (t  ,s  )  ( i  =  1,3)  do  not  depend  on  x,  but  do  de¬ 
pend  on  y.  In  matrix  form 


„<ni>  = 


4> 


(j) 


0,  .  .  .  , N-l 


T(i+i> 


of  order  :  nj+p  x  nj 


7 


s(j+i) 


j=  1 


.  .  ,N-1 


s(j+1) 


of  order  nj+]_  x  nj-l 


Suppose  we  found  and  such  that 


<P 


<i>  -  b(5>  <,,<0>  *  a<5>  A 


where 


B(^  :  n.  x  (2N+1) 

3 

A(j)  :  n.  x  (2N+1) 

3 

B(D  =  T(1)  A(1)  =  S(1) 


B(0)  =  I  A(0)  =  0 


1) 


then : 
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The  coefficients  of  the  scheme  are  c btained  Ly  the  demand  that 
it  is  accurate  for  3-polynomial  lead.  Since  the  7-point  formula  is 
accurate  for  8  polynomials  the  order  of  accuracy  is  not  changed. 
The  9-point  scheme  can  be  expressed  as: 


where : 


^(j+3/2)  =  “(j+3/2) 

~(j+3/2)  , 
"(j+3/2)  . 


4,<j>  +  s(3  +  3/2) 

n j  +3/2  x  nj 
n j  +  3 / 2  X  nj-l 


<P 


(j-1) 


and  1)  can  substituted  from  the  equation  obtained  by 

the  7-point  formula: 

„<3>  =  b(3>  .  A(j>  v'-1* 

,ip<j*3/2)=i(3*3/2)(B(j)(()(0)tA(j)v(-l))t-(3*3/2)(B(j-1)v(0) 

♦a(3-1)*(-1)) 

=  B<i»3/2)(t(°)M(j*3/2)#<-l) 


where : 
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g( j+3/2)  =  ^(j+3/2)g(j)  +  g(j+3/2)g(j-l ) 
A<j+3/2)  =  “(j+3/2)A(j)  +  g(j+3/2)A(j-l) 


We  now  have  2N-1  equations  for  the  2N-1  unknowns  on  the  para¬ 
bolae  line: 


g.]=4)(^)(l)=B(j)(l)tp(0)+A(j)(l)(p(_1)  j  =  1 ,  .  .  .  ,N 

gj+^^(j+1/2)ci)  =  i(j+1/2)(i)tp(0)+Aj+1/2(iK(-1)  j  =  i . 

where  B(j)(l),  A(^(l),  B(^+1/2)(l),  A(  j +  1/ 2 -1  ( 1 ) 

are  the  first  rows  of  A^  ,  B^  +  1/2),  A^  +  1/2)  respectively. 

The  projection  on  the  parabolic  line  can  be  written  in  matrix  form 


The  projection  matrices  A  and  B,  of  order  ( 2N+1 )x( 2N  +  1 )  are: 
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10  . 0 

'0  . o 

b(1)  (1) 

A(1)  (1) 

g  (  3  /  2 ) 

AC3/2)(1) 

B(2)  (1) 

A(2) (1) 

B(5/2) (1) 

^  - 

r — 1 

W 

CN» 

N 

in 

r< 

B(N-1) (1) 

* 

/i  Cl) 

B(N"^  (1) 

:(i;-i/2)(D 

/■"s 

i — 1 

CD 

A(:<  }  Cl) 

0 .  0  1 

—  — — i 

0 . 0 

__ 
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4 .  Solution  of  the  elliptic  problem 

In  the  elliptic  region  we  .se  ei'he:  a  5 -point  (I)  or  a  9-point 
difference  scheme. 


Ay"  ip  ,  «  tp(x,v+  '  v)  -  .  •: 1  , v )  *  ip(x,y-£y) 
j  y 


,  n  -  «:4>.  •  +  4>- 

1 , 1  +  1  ’  -  -  > 


Ax  4>xx  «  43(x+Ax,y)  -  24>(x,y)  +  4)(x-Ax,y) 


=  IL  .  ,  .  -  2w .  .  +  w . 

V1+1,J  V1,D 


The  5-pcint  difference  scheme  for  the  Triccmi  eouaticn  yip  -4>  : 

xx  yy 


for  the  point  (i,j)  is: 


where : 


-  2cp  +  ip  )  -  (ip  .  -  2ip  +  4>;  •_-) 

1  i  J  1  A  1  j  1  5  J  1  -*-  f  .  ^  5  _  ^ 


y  =  jAy  j  =  -M  +  1,  -M  +  2 , .  .  .  ,  1 


Ihia  scheme  is  accurate  for:  l,x,y,xy,3x"  +  y 


1)  We  build  a  3-point  difference  scheme  by  taking  the  lire 


polynomials  which  solve  Tricomi's  equation,  i.e.: 
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2  <'<3  2  4  3  4 

l,x,y,xy,lx  +  y" ,  '  +xy  ,  6x  y  +  y  ,2x  y  +  xy 


The  demand  that  the  difference  scheme  based  on  the  points: 

Cx.  ±  Ax,y.  i  Ay)  ,  (x.  +  Ax,y.)  ,  (x.,y.  +  Ay)  and  (x.,y.) 

13  131]  13 

accurate  for  those  functions  leads  to  a  dependent  system  of 

equations.  Therefore  one  polynomial  has  to  be  exchanged  with  or.e 

with  a  higher  degree. 

An  independent  set  of  functions  are: 


2  3  3  j  ?  14  3  14  U  '  3  0 

l»yjxy»3x  +  y  ,x  +xy  ,6x“y  +  y  ,2x  y  +  xy  ,15x  +  33x“y  +  2y 


The  scheme  is  accurate  for  P  =  as  well,  hence  it  is  accurate  fo 
9  functions. 

Both  types  of  the  difference  schemes  give  a  system  of  equations 
which  is  triangular  in  blocks: 


14 


with  all  the  blocks  having  the  same  size:  (2N+1)  x  (2N+1) 


^  -M+l  j  j  <  1,  are  matrices  of  either  the 

5-point  or  9-point  difference  schemes  including  identity  equa¬ 
tions  for  the  boundary  points  (-l,jAy),  (+l,jAy)A^^  and 

(0)  (-M) 

B  are  the  projection  matrices  (section  3),  f  is  the 

boundary  condition  on  y  =  -1,  and  f^^Cl)  and  f^^(2N+l)  are 

the  boundary  conditions  on  x  =  -1  and  x  =  +1,  accordingly. 


f(j)(l)  =  ip(-l,  j  Ay) 

f(j)(2N  +  l)  =  «p(+l,jAy))  -M  +  1  S  j  S  1 
fj(i)  =  0,  2  *  i  f  2N 

. 

-(0)  ,  .  .. 
f  =  g  -  vector  of  solution  on  F^. 


Numerical  results 


-  15,  - 


The  mixed  problem  is  solved  in  3  steps: 

a)  Computation  of  the  projection  matrices  and 

b)  Solution  of  the  elliptic  problem  with  A^^ip^  ^  +  B^^ip^^g 

as  boundary  condition  on  y  -  0 . 

c)  Solution  of  the  hyperbolic  problem  using  the  results  of  if/ ^ 

( -1 ) 

and  tp  from  (b),  ana  proceeding  with  the  7-point  formuia: 


To  see  the  influence  of  the  projection  condition  on  y  =  0  the  ellip¬ 
tic  problem  is  also  solved  with  exact  boundary  conditions  on  y-0, 
i.e.  A^^  =  0,  =  I  and  f^^  =  ip(x,0)  is  substituted  in  the 

system  (d);  the  hyperbolic  problem  is  also  solved  with  exact  val- 
ues  for  ip  and  ip  .  (e) 

The  folloiwng  3  examples  are  considered: 


i)  ip  (x,y)  =  15x4  +  30x^y^  +  2y 


Note:  For  this  function  the  9-point  formula  in  the  elliptic 

region  is  exact. 

4  2  4  7 

ii)  ip0(x,y)  =  -?lx  y  -  2 lx  y  -  y 


n+i 


iii)  ipQ(x,y)  -  coshx ( y  +  £ 


(  3n  + 1 )  3n  (  3n-  2  )  (  3n-  3  )  .  .  .4.3 

n- 1 
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Conclusions 

The  order  of  convergence  of  the  5-point  formula  alone  is 

2 

0(Ax  )  and  the  order  of  convergence  of  the  9-point  formula  alone 
4 

is  0(Ax  )  (see  results  (d),  table  1-3).  The  results  of  (e)  show 

3 

that  the  order  of  convergence  m  the  hyperbolic  region  is  0(Ax  ), 
approximately : 

0(Ax3  +  2/3)  «*  0(Ax3*Ayj)  (j>0) 

(estimate  of  error,  see  [3]). 

The  schemes  used  in  the  hyperbolic  region  are  more  accurate  than 

the  5-point  formula  and  less  accurate  than  the  9-point  formula, 

used  in  the  elliptic  region.  Therefore,  when  the  mixed  problem  is 

solved  with  the  projection,  the  maximal  absolute  error  is  obtained 

in  y  <  C  in  case  of  the  5-point  scheme  and  the  order  of  con- 
2 

vergence  is  0(Ax  )  in  the  whole  region,  and  the  maximal  absolute 
error  is  achieved  in  y>0  when  we  use  the  9-pcint  scheme  in  the 
elliptic  region  and  the  order  of  convergence  is  0 ( Ax ° Ay . )  (j>0) 
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